Variance in heat tolerance in bumble bees correlates with species geographic range and is associated with several environmental and biological factors

Abstract Globally, insects have been impacted by climate change, with bumble bees in particular showing range shifts and declining species diversity with global warming. This suggests heat tolerance is a likely factor limiting the distribution and success of these bees. Studies have shown high intraspecific variance in bumble bee thermal tolerance, suggesting biological and environmental factors may be impacting heat resilience. Understanding these factors is important for assessing vulnerability and finding environmental solutions to mitigate effects of climate change. In this study, we assess whether geographic range variation in bumble bees in the eastern United States is associated with heat tolerance and further dissect which other biological and environmental factors explain variation in heat sensitivity in these bees. We examine heat tolerance by caste, sex, and rearing condition (wild/lab) across six eastern US bumble bee species, and assess the role of age, reproductive status, body size, and interactive effects of humidity and temperature on thermal tolerance in Bombus impatiens. We found marked differences in heat tolerance by species that correlate with each species' latitudinal range, habitat, and climatic niche, and we found significant variation in thermal sensitivity by caste and sex. Queens had considerably lower heat tolerance than workers and males, with greater tolerance when queens would first be leaving their natal nest, and lower tolerance after ovary activation. Wild bees tended to have higher heat tolerance than lab reared bees, and body size was associated with heat tolerance only in wild‐caught foragers. Humidity showed a strong interaction with heat effects, pointing to the need to regulate relative humidity in thermal assays and consider its role in nature. Altogether, we found most tested biological conditions impact thermal tolerance and highlight the stages of these bees that will be most sensitive to future climate change.


| INTRODUC TI ON
Across much of the northern hemisphere and South America, where the ~270 bumble bee (Hymenoptera: Apidae: Bombus spp.) species are native, numerous flowering agricultural crops and wild flowering species have evolved to be heavily dependent on bumble bees for their pollination (Goulson, 2010).These services, however, are under threat given documented declines in these bees over the last century (Cameron et al., 2011;Goulson et al., 2008;Soroye et al., 2020).
In Europe and North America, 26% of bumble bee species now are listed as threatened or endangered (Cameron & Sadd, 2020), and bumble bee species richness has declined 30% in the Northeastern United States since the late 1800s (Bartomeus et al., 2018).While there are many stressors impacting bee populations (e.g., pathogens, habitat loss, and pesticides), climate change is a major factor linked to range shifts and biodiversity loss in bees and other pollinators (Soroye et al., 2020;Vasiliev & Greenwood, 2021), with several studies in bumble bees exemplifying these effects.
Bumble bees are particularly vulnerable to increasing temperatures associated with climate change.Hamblin et al. (2017) found bumble bees to be most susceptible to heat among 15 tested bee species in North Carolina, U.S.A., and Pardee et al. (2022) found that bumble bees were less abundant during warmer years compared to 28 other bee genera in a montane region in Colorado, U.S.A.Although bumble bees could benefit from some warming as it might reduce their expenditures toward maintaining endothermy (Jackson et al., 2020), bumble bees have been found to decrease foraging and their colonies to have decreased growth with high heat (Gérard et al., 2022;Hemberger et al., 2023;Vanderplanck et al., 2019).Effects of climate change in bumble bees have been well documented using large-scale geographic and climate modeling approaches.In Europe and North America, climate has been found to be a better predictor than habitat in explaining declines, and regions with greater increase in temperature have experienced greater declines in bumble bee diversity (Soroye et al., 2020;al-though see Guzman et al., 2021).Bumble bee species have also experienced range contraction in response to climate warming, with southern ranges receding northward without a corresponding shift at the northern edges of their ranges (Kerr et al., 2015).Bumble bees may also be shifting altitudinally with climate, potentially creating isolated populations that are vulnerable to extirpation (Kerr et al., 2015;Pyke et al., 2016).Forecasting models taking into account climatic niches and dispersal abilities of bumble bees predict that many species will fail to disperse into new climatically suitable conditions with projected climate change, leading to declining populations (Sirois-Delisle & Kerr, 2018).Not all bumble bee species have declined equally in the face of recent climate and land use changes (Cameron et al., 2011;Colla et al., 2012;Grixti et al., 2009).Determining the factors that make species resilient requires a better understanding of the factors that limit species distributions.Bumble bees are cold-adapted with the highest species diversity in cool, temperate-boreal, and montane habitats, although a few species are endemic to deserts and tropical lowlands (Hines, 2008;Williams et al., 2014).Recent studies point to cold tolerance being a limiting factor on altitudinal and geographic distribution, and cold tolerance is also likely under selection (Jackson et al., 2020;Keaveny et al., 2019Keaveny et al., , 2023;;Maebe et al., 2021;Oyen et al., 2016;Pimsler et al., 2020).However, the global distributions of bumble bee species appear to be further influenced by species-level variation in heat tolerance (Table 1): Martinet, Dellicour, et al. (2021) performed an analysis of thermal tolerance on males of 39 bumble bee species across Europe and North America and found interspecific variation in heat stress tolerance that matches species distributions.Arctic and boreal species had lower resistance to hyperthermic stress than temperate and Mediterranean species.The high heat tolerance in the most southern distributed Mediterranean species Bombus terrestris (Linnaeus) is particularly striking given that this species is generally expanding its range.Oyen and Dillon (2018) likewise found that the lower altitude Rocky Mountain species Bombus huntii Greene had higher heat tolerance than Bombus sylvicola Kirby, which inhabits higher altitudes.Gonzalez et al. (2022) applied thermal tolerance data from other studies and noted a relationship between heat tolerance and several climatic variables in five North American bumble bee species.Differential heat tolerance in bumble bees may be conferred through differences in molecular heat shock response (Blasco-Lavilla et al., 2021;Kuo et al., 2023;Pimsler et al., 2020), interspecific differences in heat shunting by means of counter-current exchange between the thorax and abdomen (cf.Heinrich, 1976;Heinrich & Vogt, 1993), or through other physiological adaptations such as shifting metabolism or potential for evaporative cooling (Johnson et al., 2023).
The susceptibility of bumble bees to climate is likely to vary not only across species but by life history and condition within species.Prior studies examining species heat tolerance have shown high variance across individuals within species (Martinet, Dellicour, et al., 2021;Oyen & Dillon, 2018), suggesting that other biological factors interplay with thermal physiology to impact tolerance.Bumble bees, as social species, exhibit variation across individuals in physiology and behavior depending on caste and stage of their life cycle (Amsalem, Grozinger, et al., 2015).Queen bumble bees, for example, are known to shift their nutritional reserves in preparation for diapause and upon ovary activation (Amsalem, Galbraith, et al., 2015;Treanore & Amsalem, 2022), and workers in late or queen-less colonies can also lay eggs (Goulson, 2010).
There may also be selection for bumble bees in different roles to have different heat tolerances based on differences in thermal

Ecophysiology
TA B L E 1 Summary of the literature on bumble bee heat tolerance for the parameters examined here.Provided are the parameters studied, whether or not there was an effect ("Eff?") on heat tolerance, the bumble bee species studied, how the parameters relate to heat tolerance, and the citation of the study.

Species
Yes 39 spp.From three continents Higher THS in species from warmer habitats exposure.Bumble bees typically nest underground where they tend to be buffered from the elements, although some species tend toward surface nesting (Liczner & Colla, 2020).Males are exposed to environmental conditions more than workers, as workers can escape the heat in nests (Heinrich, 2004).Among workers, foragers are more exposed to heat stress than nurses, which usually reside only in the nest.Queens in general receive little exposure to heat, as they are typically produced in late summer where they initially feed in natal nests, subsequently receiving the most exposure when they leave the nest to mate (~1-2 weeks postemergence; Röseler & van Honk, 1990), after which they dig themselves a hibernaculum underground for overwintering.They will emerge again in the cool of spring to forage for provisions for their nest, thereafter remaining protected in underground colonies (Gardner et al., 2007).
There have been a few studies examining how thermal heat tolerance in bumble bees varies with sex and caste (Table 1).Age-related variation in thermal tolerance has been documented in a number of holometabolous insects with many species showing a decline in tolerance with age (Bowler & Terblanche, 2008).Bumble bee workers tend to live from 2 to 6 weeks and can switch roles from in-nest and nursing tasks to foraging as they age (Tobback et al., 2011).A study examining the critical thermal limits of Bombus impatiens Cresson at 3, 4, and 7 days, found 3-day-old and 7-dayold workers are more heat resistant than 4-day-old workers (Oyen & Dillon, 2018).The study, however, did not test the effect of older ages (>1 week) on thermal tolerance.
Bumble bees are facultative endotherms and can maintain constant body temperatures through a wide range of air temperatures (Church, 1960;Goulson, 2010;Heinrich, 2004).Nevertheless, by Bergmann's rule, the small surface to volume (S/V) ratio of larger bees should reduce their thermoregulatory capabilities for dissipating heat (Heinrich, 2004), while increasing their ability to retain heat in cold conditions.Given evolutionary tradeoffs in body size in insects, however, Bergmann's rule may not always apply (Lozier et al., 2021;Shelomi, 2012).Prior work found that CT max increased with mass in reared workers of B. huntii (Oyen et al., 2016), but not in wild workers of B. bifarius, B. sylvicola (Oyen et al., 2016), B. impatiens (Oyen & Dillon, 2018), or B. terrestris (Maebe et al., 2021).Martinet, Dellicour, et al. (2021) found only weak explanatory power of body mass on time to heat stupor at 40°C among males of 39 bumble bee species.Recent studies have found a decrease in body size in workers in the last 100+ years in certain bumble bees (Nooten & Rehan, 2020) and other bee populations (Oliveira et al., 2016).Although developmental plasticity could be involved in this, it is possible that climate change may be implicated via natural selection if larger bees are more negatively impacted by heat.
More evidence is needed to understand if body size indeed plays a role in heat stress tolerance and long-term shifts in body size in bumble bees.
Nutrition may also impact thermal tolerance (Vanderplanck et al., 2019) and, if so, this can both explain natural variance in heat tolerance and present a solution to mediate effects of thermal stress.There have been mixed results on the effects of starvation on thermal tolerance with some studies finding a role (Blasco-Lavilla et al., 2021;Quinlan et al., 2023) and another seeing no correlation (Oyen & Dillon, 2018).
In we examine the impact of humidity levels on resulting heat stress tolerance.Humidity and temperature interactions may be important as high humidity can reduce the potential for evaporative cooling in insects as it reduces the hygric gradient between the bee and environment (Church, 1960;Prange, 1996).While evaporative cooling is thought to not be a major factor in most bees (Heinrich, 2004;Johnson et al., 2023), in bumble bees high humidity has been found to reduce cooling during flight (Church, 1960).
Altogether, our data reveal that many biological factors impact thermal tolerance in bumble bees, providing a baseline for understanding points of vulnerability and for guiding solutions to better manage bumble bees under climate change.

| Bee specimen sources
Bees for our thermal trials were sourced either directly from the wild (collected as foraging queens, workers, or males), drawn from laboratory-reared colonies founded by wild-caught queens, or from commercially purchased colonies (thus from laboratory-reared queens).While we sought to balance replicates in the design as much as possible with our sampling effort, sample size of field-collected bees was somewhat uneven across groups due to limited field availability.A total of 379 bees were tested for these assays with sample sizes and species information contained in Table A1, with a more detailed table in Scholarsphere (https:// doi.org/ 10. 26207/ gqy8-e215), and apparent in Figure 1A.

| Wild bees
A bee is defined as wild if thermal tolerance was tested within 24 h of its collection from nature.sucrose and 50% invert sugar with amino acid supplement (Amino-B Booster; Honey-B-Healthy Inc., Cumberland, MD) to provide bees with essential amino acids.

| Laboratory-reared bees from commercial colonies
For some of our experiments, we assayed workers, males, and newly emerged queens (we use the term "queen" here to include both mated or unmated and those that have or have not started a nest) produced from commercial, research-grade B. impatiens colonies (sourced from Biobest, Canton, MI).Commercial colonies were kept under the same conditions as laboratory-reared wild queen colonies.

| Time to heat stupor protocol
We used time to heat stupor (THS) as a measure of thermal stress across all bees (Martinet, Dellicour, et al., 2021).Symptoms of heat stupor include a loss of neuromuscular function, which is characterized by uncoordinated movements, the inability right itself when flipped on its back, twitching in the extremities (e.g., tarsi), and heat coma (cf.Martinet, Dellicour, et al., 2021).We used THS instead of CT max (critical maximum temperature, the temperature at which bees reach heat stupor following an incremental temperature ramping rate from ambient temperature; Oyen & Dillon, 2018), as we sought to test thermal endurance without the confounding factor of ramping rate.We felt THS would be more akin to the temporal response upon entering a hot environment from a cool nest.Moreover, THS has been useful for discriminating bumble bee species-specific tolerances in the past (Martinet, Dellicour, et al., 2021;Zambra et al., 2020).Static (THS) and dynamic (CT max ) assays of heat tolerance have been found to be comparable and thus reliable methods of assessing heat tolerance in Drosophila (Jørgensen et al., 2019) and thus these methods should yield similar results in bumble bees, however, we recognize that the method applied will likely affect magnitudes and variance of response.Trials were limited to 3 h to avoid confounding effects of starvation on THS.We thus consider bees reaching 3 h as "heat resistant."In our preliminary THS trials, we used 40°C as our critical temperature.This temperature has been used in previous studies to induce heat stupor across bumble bee species (Martinet, Dellicour, et al., 2021;Zambra et al., 2020) and is above the threshold temperature commonly used to define a heat wave (Xu et al., 2016).However, most bees in our preliminary study were reaching 3 h without entering heat stupor, so we decided to use a higher temperature of 43°C that would capture the full range of responses to heat stress within the trial period.
We used a water bath (Benchmark, MyBath 4L, model H2004) for these trials, to retain consistent heat and humidity (constant at 90% RH) across periodic (5 min) THS checks.Bees were taken from their holding container/colony and transferred to preheated (43°C ± 0.5°C) cylindrical glass vials (volume = 35.35cm 3 ) (Genesee Scientific, Flystuff, 32-117BF) placed in a rack in a water bath (Benchmark, MyBath 4L, model H2004) (90% RH).Glass vials were weighed down with hex nuts (29.5 g) glued to the bottom of the vials and were cotton-stoppered to allow for gas exchange and prevent CO 2 buildup from respiration.Bees were checked every 5 min for symptoms of heat stupor by briefly lifting them from the water bath and flipping or shaking the vial as needed to assess if they could remain upright.

| Cross-species analysis
We performed thermal tolerance assays on collected workers, (A) Interspecific differences in THS across all variables, highlighting differences by species, sex/caste, and between wild (W) versus laboratory-reared (L) bees for each caste-species combination.For ease of viewing, significance letters in 1a represent pairwise differences among species/origin (Tukey HSD) within each caste.(B) Interspecific differences among workers (wild and laboratory-reared bees pooled) highlighting significance values and how THS of species directly relates to their latitudinal averages of GBIF records, mean temperature (°C) across this range for June, July, and August (indicated with Xs), and habitat preferences across a forested hill to open and developed valley landscape in Pennsylvania (Gratton et al., 2023).Significance letters on the left of the slash represent pairwise differences (Tukey HSD) between workers of different species and letters on the right of the slash represent pairwise differences (Tukey HSD) between species from the full model, which includes workers, males, and queens in the analysis.(C) Sex/caste differences in laboratory-reared Bombus impatiens, which had the best sample sizes to observe caste differences, and their significance.Sample sizes for each condition are indicated in parentheses.bimac, Bombus bimaculatus; gris, Bombus griseocollis; imp, Bombus impatiens; perp, Bombus perplexus; v/s, Bombus vagans/sandersoni.
(Table A1).B. vagans and B. sandersoni are less common, are hard to tell apart morphologically, and occupy the same ecological niche in Pennsylvania higher elevation forests (Gratton et al., 2023) For these two models (additive and interactive), we also calculated Akaike's information criterion, corrected for small sample size (AICc), using the AICcmodavg package (Mazerolle, 2017).While these models reveal the main effects, we also conducted separate analyses for visualization purposes on subsets of the data.In particular, we pooled wild and laboratory-reared workers within each species and compared THS among species.We also compared laboratory-reared B. impatiens workers, males, and queens, for which we had the most individuals of each sex/caste.We chose to exclude wild bees from this analysis due to potential confounding effects of origin on THS.
To assess species niche (i.e., latitude, climate, habitat type) relative to heat tolerance results, we used historical, georeferenced records of each Bombus species from GBIF (GBIF.org, 2023a(GBIF.org, , 2023b(GBIF.org, , 2023c(GBIF.org, , 2023d, 2023e, 2023f), 2023e, 2023f).For latitude, we filtered specimen records from the United States and Canada and only used specimen occurrence records east of 100.00° longitude to focus on eastern North America.From this, the 1% tails of the distribution were excluded to focus on the core range of each species.Within each of these ranges, we obtained average summer temperatures for the months when workers are primarily foraging-June, July, and August-aggregated over 30 years .Climate data was extracted from PRISM Climate Group (2020) using the prism package (Hart & Bell, 2015).We used separate ANOVA models followed by Tukey's HSD tests to assess both species differences in mean latitudinal distributions and summer mean temperatures.As PRISM data does not include climatic data for eastern Canada, Canadian specimens were not included in the temperature analysis.We also compare our results to data on altitudinal and habitat niche (forested hills, open valleys, and the transitional habitat between these) of each of these species in Pennsylvania (Gratton et al., 2023).

| Colony
We tested the effect of colony origin on THS in workers from five commercial colonies (n = 3-5 bees/colony) and three colonies reared from wild queens (n = 13-16 bees/colony; Table A1), using a two-way ANOVA followed by a Tukey's HSD test, testing the effect of both colony and wild versus commercial origins.There was no significant variation in thermal tolerance with colony and thus we did not include colony identity in future models.

| Age
We compared THS among B. impatiens queens of different ages, drawn from three different commercial colonies: 4-day-old (n = 17), 7-day-old (n = 4), 10-12-day-old (n = 9), and 12-32-day-old (n = 17) (Table A2).For the first three age groups, callow queens were removed from their natal colonies daily, and individuals from the same colony were kept in a container together and provisioned with artificial nectar and pollen ad libitum until they reached the desired age.The 12-32-day old cohort was removed from natal colonies at 2-20 days after adult eclosion and kept together prior to sampling.
We also compared THS among two different age groups of B. impatiens workers.Seven-day-old (n = 12) and 14-day old workers (n = 13) were sampled from three laboratory colonies reared from wild-caught queens.Callow workers were removed from their natal colonies and placed in mini-colonies (a few individuals together per container) provisioned with artificial nectar and pollen ad libitum until they reached the desired age.The effect of age on mean THS for both queens and workers were tested separately using ANOVAs followed by Tukey's HSD tests.

| Ovarian activation in queens
We sought to test whether ovarian and physiological state may im- were used as controls.We divided treated and untreated bees into two subgroups in a fully crossed design: 2-day post-CO 2 narcosis/2day control (n = 8 each) and 7-day post-CO 2 narcosis/7-day control (n = 9 each) (Table A2).After queens went through THS trials, queens were dissected, and ovary activation was measured by averaging the length of the largest terminal oocytes from each ovary, plus the largest remaining oocyte from one side (cf.Amsalem & Grozinger, 2017).
We used ANOVAs followed by Tukey HSD tests to assess whether there are differences in mean THS between CO 2 narcosis treatments and ovary activation (average terminal oocyte length).To further test the effects of ovarian activation on heat stress response, we performed a linear regression model comparing THS against average terminal oocyte length across all samples.Finally, some queens had started laying eggs, so we performed ANOVAs to test the effects of presence/absence of brood on mean THS and examined whether bees with brood indeed had more ovary activation.

| Body size
To test the effect of body size on heat tolerance, intertegular distances (ITD) were measured in 71 B. impatiens workers that underwent THS trials.ITD is a commonly used metric to assess body size in bumble bees and other bees and is considered a reliable metric of body size (Cane, 1987;Hagen & Dupont, 2013).ITD was measured with a dissection microscope ocular reticle and converted to millimeters for analysis.We did separate analyses for workers reared from wild-caught queens (n = 55) and for wild workers (n = 17).For each group of workers (laboratory-reared and wild), we used a simple linear regression model to test the effect of ITD on THS.The data for laboratory-reared workers appeared somewhat curvilinear, so we also ran a second-order polynomial regression on this data.We then compared this model to our linear regression for laboratory-reared workers using AICc.We retained the linear regression model since it had the lower AICc (172.7 compared to 176.0).

| Humidity and temperature interactions in THS
To test the interaction between temperature and humidity, we first conducted an experiment on commercial B. impatiens workers (n = 4 colonies; three replicates per colony/condition; n = 12 bees/condition) at three different humidity conditions: 20%, 50%, and 70% RH (range ± 4% during trials).These humidities were chosen to span dry, average, and higher humidity conditions.Bees were taken from their natal colony, placed into separate dry incubators (Thermo Scientific, 371 L) set at 43°C (±2°C during trials), and checked every 5 min for heat stupor.
As the temperature fluctuated quite a bit during the first experiment, likely due to our regular opening of the door, we ran a second experiment using a single device for all specimens, seeking to control environmental conditions more.Given the results from the first experiment, we focused only on the highest and lowest humidities (20% and 70%).For this, we drew 10 workers randomly for each condition (20% [17%-29%] and 70% [67%-77%]) from a single commercial parent colony that was generating only workers, and we also drew 10 queens (aged 10-to 21-days-old) from a different commercial colony that was in a queen-producing stage.Bees were moved into a small (23.5 cm × 27.5 cm × 36 cm) incubator (Vevor, model: XHC-25) that was set to 43°C (actual temperature: 42.9-43.5°Cat 20% and 42.4-42.7°Cat 70%), set within a walk-in incubator at similar humidity to the treatment and at 28°C, and assessed for THS.
We compare these data with THS results from B. impatiens workers (n = 69) run at a stable 43°C and 90% RH in a water bath.Note that only for humidity trials did we use incubators; water baths were used for all other trials.
For the humidity trial analysis, given that much of our data were truncated to 3 h due to the great number of heat resistant bees at low humidity, we ran a GLM binomial regression to assess the proportion of bees that entered heat stupor as a function of humidity (20% vs. 70% RH) and caste (workers and queens).We also ran an ANOVA followed by Tukey tests on both humidity datasets separately, testing the effects of humidity treatment for the first experiment and both humidity and caste for the second.

| RE SULTS
3.1 | Cross-species analysis (species, caste, and origin) of heat tolerance

| Species
We found a strong effect of species (p < .001) on mean THS in our full model (  1A).

| Origin of bees
We found a significant difference between laboratory-reared and wild bees overall (p < .001),an interactive effect of species and origin (wild-vs.laboratory-reared) (p = .002)(Table 2), and a trend for interactions of origin with sex/caste and species (p = .062).Across species, wild-caught workers had higher THS on average than laboratory-reared workers with the exception of B. griseocollis, where the opposite pattern was observed (Figure 1A).

| Body size
We found larger body size was associated with decreased THS for wild B. impatiens workers (R-squared = .26,F 1,15 = 5.29, p = .036),but there was no effect of body size on laboratory-reared workers (Figure 4).On average, wild workers had a larger body size than laboratory-reared workers (p = <.001).

| DISCUSS ION
Prolonged hyperthermic stress in bumble bees can have lethal as well as sublethal effects on bee health, such as limiting forage time and flight ability (Hemberger et al., 2023), and impacting colony growth and reproduction (Vanderplanck et al., 2019).The frequency and duration of heatwaves across the globe are predicted to increase due to climate change (Domeisen et al., 2023), and thus many bumble bee populations are more likely to experience hyperthermic stress in the future.To better understand how these bees will respond to this change, we have provided data toward understanding how different species respond to heat stress and which physiological, morphological, and life-history factors make them more vulnerable.We found that species, caste and sex, origin, age, ovary activation, body size, and relative humidity all play a role in how the bumble bee species in our study region respond to heat stress, but to different extents.
Most notably, we found that bumble bee species vary in heat tolerance in a way that matches their habitat distributions, suggesting that heat sensitivity likely impacts their range.We find that queens are most vulnerable, but that their vulnerability is dependent on age and reproductive status.Our field to laboratory comparisons suggest that some of the exceptional variability we find in heat tolerance in the field may be due to physiological status or environmental conditions.Our data further highlight the strong effect humidity plays in thermal resilience and the need to consider humidity and temperature interactions in future work.Our results identify points of vulnerability in these bees to consider when addressing decisions and policy regarding bumble bee conservation, especially in the context of climate change.Martinet et al. (2015) and Martinet, Dellicour, et al. (2021), in comparing bumble bees across the globe, found differences in heat tolerance that reflected their respective ecoregions.Other work in montane regions have found that bumble bee altitudinal ranges are correlated more with cold tolerance than with heat tolerance (Pimsler et al., 2020).

| Thermal tolerance of species tracks their distribution
Large-scale climate modeling lends support to heat vulnerability driving shifting ranges, range contraction, and species loss, suggesting natural variation in heat tolerance may be a reason that some species are particularly vulnerable (Kerr et al., 2015;Sirois-Delisle & Kerr, 2018;Soroye et al., 2020).We find a clear trend in inferred heat tolerance that matches the temperature of zones these species occupy.Species with higher mean THS have more southerly median latitudinal distributions, warmer climatic distributions, and prefer valleys to forested and higher elevation habitats.Specifically, we found that B. griseocol- Bumble bees may also vary in nesting position relative to the ground surface, which can impact the ability to escape heat and the overall exposure to ambient temperatures.For the species studied here, B. impatiens is a primarily underground nester, B. bimaculatus, B. vagans, and B. perplexus are thought to be mostly underground nesters with some surface nesting, and B. griseocollis is primarily a surface nester (Husband, 1977;Plath, 1922Plath, , 1927Plath, , 1934)); thus, B. griseocollis would have nests that are less buffered from heat waves.B. griseocollis did show slightly higher heat tolerance than B. impatiens, which was not well matched by differences in the thermal occupancy between them.

F I G U R E 2
Differences in mean time to heat stupor relative to age in Bombus impatiens queens and workers.Individuals are age-staged from commercially reared B. impatiens colonies.Significance letters represent differences within each caste (queens or workers).Numbers in parentheses are sample sizes.min, minutes.
These results suggest that some species in a community are more vulnerable to heat stress than others, and that the increased frequency and duration of heat events under climate change may disproportionately impact bumble bee species with cooler climatic distributions that typically inhabit forests and more northerly latitudes.For example, the relatively high heat tolerance of B. impatiens and B. griseocollis may be one reason these species' populations are increasing proportionally in the Eastern United States and Midwest in the last several decades (Averill et al., 2021;Jacobson et al., 2018;Novotny et al., 2021) compared to the species that occupy cooler areas.B. vagans and B. sandersoni may be especially impacted since they not only have cooler climatic distributions, but have suffered recent population declines in parts of their ranges (Colla et al., 2012;Grixti et al., 2009;Jacobson et al., 2018).These results can help inform conservation management decisions by prioritizing species that are particularly sensitive to extreme heat events and are already of conservation concern.

| The environment influences variance in thermal tolerance
We found that wild bees have higher tolerance to hyperthermic stress than bees reared in the laboratory (with the exception of B. griseocollis), indicating that environmental conditions may be involved in heat stress resilience.In contrast, Gonzalez et al. (2022) found laboratory-reared bumble bees to have slightly higher thermal tolerance.Although it would require further testing, laboratory-reared bees were fed honey bee-collected pollen which did not allow the bees to seek out their nutritional optima.For example, honey bee pollen is lower in protein:lipid ratios compared to what bumble bees obtain when foraging in the wild (Vaudo et al., 2018).This could have reduced food quality and have negative effects on heat resilience at the colony (Vanderplanck et al., 2019) and individual (Quinlan et al., 2023) level.Other factors besides nutrition may also vary between field and laboratory, such as pathogen levels, optimal hydration, and ability to attain optimal nest conditions: these factors need further study to fully understand their impact on thermal tolerance.Furthermore, laboratory-reared bees may have included nurse bees, while wild-collected bees are only foragers.While wild bees will have experienced greater temperature fluctuations during their lifetimes than laboratory-reared bees, recent studies demonstrate that prior experience of thermal stress does not result in increased thermal tolerance through priming (Quinlan et al., 2023;Sepúlveda & Goulson, 2023).

| Sex and caste impact thermal tolerance
Across all species, and when analyzing B. impatiens alone, queens had lower THS than workers and males.The relatively high sensitivity of queens to hyperthermic stress in the laboratory appears to make them particularly vulnerable to extreme heat compared to workers and males.This effect, however, was only observed in our water bath trials where relative humidity was 90%.At 70% RH, queens were not different than workers, indicating an interaction between caste and humidity.
In our assays, males were slightly less heat tolerant than workers, although the significance was marginal, and this showed a trend toward being dependent on species.Prior work found no difference in CT max between males and workers in three species in the Rocky Mountains (Oyen et al., 2016) and no differences in CT max between worker and males among four Bombus species native to Columbia (Gonzalez et al., 2022).There is temporal overlap between male and queen emergence, both leaving the nest in the heat of late summer, putting them at similar risk of encountering the same extreme heat events.While males were less sensitive to heat stress (had higher THS) than queens, sperm viability has been documented to decline in male bumble bees exposed to elevated temperatures that fall far below the temperatures tested here (Campion et al., 2023;Martinet, Zambra, et al., 2021), thus cross-generational effects must also be considered in their heat adaptation and vulnerability.

| Age effects on thermal tolerance depend on caste
Our B. impatiens trials show that queens differ in their thermal tolerance by age.Four-day-old queens are the most vulnerable to heat stress, followed by a steady increase in THS as they age until around 12 days where their levels reached the higher levels typically observed for workers, after which they decline to the levels we typically observed in wild queens/queens.This suggests that queens change physiologically as they age, which impacts their response to heat stress.Oyen and Dillon (2018) similarly found that 4-day-old workers were less heat tolerant than 3-or 7-day-old workers, suggesting this may be a vulnerable stage.The peak of heat tolerance in queens at 12 days is intriguing in its biological correlates: Queens are thought to leave the natal nest to mate and find a hibernaculum at 7-12 days-old (Tasei et al., 1998;Treanore et al., 2021) and then likely enter cooler substrates for diapause thereafter.Thus, thermal tolerance may align with the period when queens are most exposed to thermal stress.
We have also expanded on data by Oyen and Dillon (2018) by examining the effects of age in older worker bees.We found no effect of age between 7-and 14-day olds on heat tolerance, despite 14-day olds being near the point of mortality (only one 21-day-old worker survived and was tested, but it had a similar thermal tolerance to 7-and 14-day olds).

| Ovary activation state impacts thermal tolerance
One factor that may make queens more vulnerable to heat is their investment in reproduction.Queens increase ovary activation 1-2 weeks after diapause or CO 2 narcosis (Amsalem, Galbraith, et al., 2015).Ovary activation in queens is associated with physiological changes such as increases in juvenile hormone (JH), ecdysteroids, and body fat metabolism (Amsalem, Galbraith, et al., 2015;Bloch et al., 2000).

| Humidity impacts thermal tolerance
Prior studies on bumble bee thermal tolerance have focused on temperature without mentioning potentially confounding effects of humidity (Gonzalez et al., 2022;Hamblin et al., 2017;Maebe et al., 2021;Martinet, Dellicour, et al., 2021;Oyen et al., 2016;Oyen & Dillon, 2018;Pimsler et al., 2020).We found that humidity can have strong effects on thermal tolerance and thus that humidity needs to be carefully controlled in these thermal experiments.Heat tolerance declined with increasing humidity (20%-70%-90% RH) for both workers and queens, with the effects being most notable at the highest humidity levels.A likely possibility for this association is that as air becomes more saturated at higher humidity levels, bees are less able to use evaporative cooling to bring down their body temperature.While convective cooling is thought to be the primary means by which bumble bees cool themselves, such as through shunting of hemolymph to dissipate heat through the abdomen (Heinrich, 2004), evaporative cooling is another means by which some insects survive otherwise lethal temperatures (Prange, 1996).Evaporation occurs through tracheal ventilation and in bees cooling through mouth regurgitation of nectar is purported (Heinrich, 2004;Prange, 1996).
Evaporative cooling has been observed to be more efficient in dry air than in saturated air in a number of insects, including bumble bees (Church, 1960;Prange, 1996).Bumble bee wet body mass has also been found to be more related to heat tolerance than dry mass, suggesting hydration may be important to thermal tolerance in these bees (Martinet, Dellicour, et al., 2021).The interaction between humidity and caste, whereby the reduced thermal tolerance of queens was supported only at higher humidity, may suggest that queens rely on evaporative cooling more than workers.Given these strong effects, there is a need for research on the physiological effects of humidity on thermal tolerance.

| CON CLUS ION
Our data support heat tolerance as a factor impacting where species can occur, both in terms of geographic ranges and habitat preferences, suggesting these bees are limited in part by heat sensitivity and this may play a role in their declines.We also show that queens are more vulnerable to heat stress than workers and males and especially when they are newly emerged, have active ovaries, and are producing brood.This supports a role of physiological and developmental state in not only influencing thermal response but also points males, and queens of wild and laboratory-reared B. impatiens (n = 129 bees tested), B. bimaculatus (n = 45), B. griseocollis (n = 39), B. perplexus (n = 17), and the B. vagans/B.sandersoni (VS, n = 28) complex F I G U R E 1 Thermal tolerances (time to heat stupor, THS) of different bumble bee species, castes, and in wild versus laboratory-reared bees.
and B. impatiens were more heat tolerant and occur in the warmer latitudes/climates and in Pennsylvania valleys, while the higher latitude forest-restricted species B. vagans/sandersoni were the least heat tolerant, and B. bimaculatus and B. perplexus, which have been found in edge habitat between these zones and were intermediate in their climatic and distributional means had more intermediate heat tolerance.

F
Effect of queen physiology (diapause condition and ovary activation) on time to heat stupor (THS).(A) Effect of CO 2 narcosis treatment and age on mean THS in commercially reared Bombus impatiens queens.CO 2 narcosis-treated bees represent a physiological state akin to having gone through diapause (spring queens) while those without it represent the state of pre-diapause fall queens.2-day treatments represent an age typically too early for ovaries to have developed and 7 days should be sufficient to allow ovary activation for diapaused queens.(B) Average terminal oocyte length (a measure of degree of ovary activation) by CO 2 treatment, showing that 7-day, CO 2treated queens had greater ovarian activation.(C) THS regressed against ovary activation across all CO 2 -treated bees.CO 2 was administered to queens aged 12-21 days (adjusted R 2 = .077,p = .062).(D) THS regressed against presence/absence of brood, showing that queens with brood have lower THS than queens without brood.(E) Average terminal oocyte length regressed against presence/absence of brood, showing that presence of brood is positively associated with ovary activation.Sample sizes for each condition are indicated in parentheses.

F
Time to heat stupor regressed against body size (intertegular distance in millimeters) for wild (triangles) and laboratory-reared (circles) Bombus impatiens workers.THS was correlated with body size in wild workers (adjusted R 2 = .26,p = .0363)but not in laboratory-reared workers (adjusted R 2 = −.0188,p = .95).F I G U R E 5 Time to heat stupor (THS) relative to humidity for commercially reared Bombus impatiens workers and queens.(A) Experiment 1 with workers.(B) Experiment 2 with both workers and queens.Significance from Tukey tests is indicated.Sample sizes for each condition are indicated in parentheses.*The 90% condition involved water baths while the remainder of the treatments were performed in air incubators.Differences in results by experiment are likely due to the way temperature was calibrated for each experiment which resulted in slightly higher temperatures for Experiment 1.
Abbreviations: Age-Q, queen age trials; Age-W, worker age trials; CO 2 Trials, CO 2 narcosis, ovary activation, and presence of brood analyses for queens; Col., colony effect analysis among laboratoryreared workers; RH1, initial humidity trial with air incubators with workers; RH2-Q, second humidity trial with queens; RH2-W, second humidity trial with workers; Sz., body size analysis among wild and laboratory-reared workers.aThenatal colony for these bees is unknown.Source: BM, Bear Meadows Natural Area; CV, Circleville Park; HV, Hyner View State Park; PS, Penn State-University Park Campus; TP, Tom Tudek Memorial Park.
, thus these were pooled to improve sample sizes.We identified B. sandersoni and B. vagans after running THS trials and found 31.6% of wild sampled bees were B. sandersoni (36.4% of workers, 37.5% of males)

Time to heat stupor assays: Bombus impatiens assays of age, reproductive state, and body size
We assessed the effect of colony origin, age, reproductive state, body size, and humidity on thermal tolerance in B. impatiens.We focused on just B. impatiens because it had more robust sample sizes for these comparisons.
heat tolerance in B. impatiens queens.Carbon dioxide (CO 2 ) 2 through a hose at low pressure for 1 min, which is enough time for bees to enter CO 2 narcosis.The hole was immediately taped over, sealing in the gas where the bees remained for 30 min.They were then placed into individual ventilated containers, where they emerged from their narcosis, and were provisioned with pollen patties and artificial nectar ad libitum.Another 17 B. impatiens queens of the same age group that did not undergo CO 2 narcosis (Amsalem & Grozinger, 2017)vation/egg development is likely already completed for CO 2 -treated bees(Amsalem & Grozinger, 2017), to examine effects of ovarian activation state.We administered CO 2 to 17 B. impatiens queens reared from two commercial colonies that were between 12-and 21-days-old age by placing bees in plastic boxes that were sealed except for a small hole in the top where we administered CO

Table 2
B. vagans/sandersoni, while B. impatiens have a higher mean THS than B. vagans/sandersoni.The interactive model had a slightly lower AICc (2477) compared to the additive model (2481), however, we report statistics from both models in Table 2.A subset model run only on workers likewise found the same patterns of significance, except that B. perplexus was not significantly different from B. griseocollis (p = .161):B. griseocollis (117.9 min ± 37.2 [SD] > B. im-

Table 3 )
. This THS order tracks the mean latitude, climatic distribution, and habitat preferences for these different species based on species distributions obtained from GBIF (Figure1B).ANOVA results for significance in time to heat stupor (THS) across conditions, including statistics for the additive MLRM on THS across species, sex/caste, and origin (top) and for the interactive multivariate model (bottom).more abundant in open valleys, B. vagans/sandersoni are more abundant in hilled forests and B. bimaculatus and B. perplexus are most abundant at the interface of these habitat zones.
(Figure1B).This matches well with the recognized differences in distribution in Pennsylvania, where B. impatiens and B. griseocollis TA B L E 2 **p-value between .01 and .001.***p-value less than .001.TA B L E 3 Results from the post hoc Tukey analyses by species (all caste/sex and workers only), sex/caste (for all species and Bombus impatiens only), and queens of different ages.Note: 10 day is 10-12 day.Bold type indicates statistical significance.Abbreviations: bimac, Bombus bimaculatus; gris, Bombus griseocollis; imp, Bombus impatiens; min, minutes; perp, Bombus perplexus; vs, Bombus sandersoni + Bombus vagans.are (González-Tokman et al., 2020) to respond to temperature(González-Tokman et al., 2020)and regulate thermal tolerance stress responses in insects (e.g., in D. melanogaster; of body size to effects in some species (Table1).We found larger body size was associated with lower THS for wild B. impatiens workers but no effect of body size on laboratoryreared workers.However, since we randomly selected workers from laboratory colonies for our trials, we likely picked some bees that were nurse bees that are restricted to in-nest activities and 4.6 | Body size effects are context-dependentBody size in principle, given the surface to mass ratio and ability to dissipate heat, should relate to thermal tolerance, with smaller bees better able to handle the heat.Prior research in bumble bees have found mixed effects of body size on thermal tolerance, from no effect